home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_13_06 / williams / matmpy.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-23  |  3.8 KB  |  133 lines
  1. /* ------------ */
  2. /* matmpy.c    */
  3. /* ------------ */
  4. #include <assert.h>
  5. #include "mpydefs.h"
  6. /* ------------------------------------------------- */
  7. /* matmpy - multiply-add/subtract operations         */
  8. /* ------------------------------------------------- */
  9. /* use:                          */
  10. /*   void matmpy(a, nra,nca, b, nrb,ncb, c, m1, m2)  */
  11. /*                             */
  12. /*   multiplies matrices a[nra][nca] & b[nrb][ncb]   */
  13. /*   and stores the result in matrix c as follows:   */
  14. /*                             */
  15. /*        (' = trnsp)  dimensions       stored     */
  16. /*    m1  product (r)  of product  m2  result (c)  */
  17. /*    === ===========  ==========  ==  ==========  */
  18. /*    AB   r =  a*b     [nra][ncb]   P   c = r      */
  19. /*    ATB  r = a'*b    [nca][ncb]  MP   c = -r     */
  20. /*    ABT  r =  a*b'   [nra][nrb]  CPP  c += r     */
  21. /*    ATBT r = a'*b'     [nca][nrb]  CMP  c -= r     */
  22. /*                             */
  23. /*   Symbols shown under m1 and m2 are defined in    */
  24. /*   header file mpydefs.h                 */
  25. /*                                                   */
  26. /*   Other parameters:                               */
  27. /*                                                   */
  28. /*      a       address of multiplicand matrix       */
  29. /*      n       ranumber of rows in a                */
  30. /*      nca     number of columns in a               */
  31. /*      b       address of multiplier matrix         */
  32. /*      nrb     number of rows in b                  */
  33. /*      ncb     number of columns in b               */
  34. /*      c       address of receiving matrix          */
  35. /* ------------------------------------------------- */
  36. # if defined(__STDC__) || defined(__PROTO__)
  37. void
  38. matmpy(double *a, int nra, int nca,
  39.     double *b, int nrb, int ncb,
  40.     double *c, int m1, int m2)
  41. # else
  42. void
  43. matmpy(a, nra, nca, b, nrb, ncb, c, m1, m2)
  44. double    a[], b[], c[];
  45. int    nra, nca, nrb, ncb, m1, m2;
  46. # endif
  47. {
  48.     double  r;
  49.     int     i, ia, ib, ic, ij, ik, incra, incrb,
  50.         j, k, mcb, mra, mrb, nfwaa, nfwab;
  51.  
  52.     /* ------------------------- */
  53.     /* Set Controls for matrix a */
  54.     /* ------------------------- */
  55.     if ((m1 == AB) || (m1 == ABT))
  56.     {                /* set for regular a    */
  57.     mra   = nra;            /* number of rows in a    */
  58.     nfwaa = nca;        /* column separation    */
  59.     incra = 1;        /* column elem incrmnt    */
  60.     }
  61.     else                        /* m1 = ATB or ATBT    */
  62.     {                /* set for a-transpose    */
  63.     mra   = nca;        /* number of rows in a  */
  64.     nfwaa = 1;        /* column separation    */
  65.     incra = nca;        /* column elem incrmnt  */
  66.     }
  67.  
  68.     /* ------------------------- */
  69.     /* Set Controls for matrix b */
  70.     /* ------------------------- */
  71.     if (m1 == AB || m1 == ATB)
  72.     {                /* set for regular b    */
  73.     mcb   = ncb;            /* number columns in b    */
  74.     mrb   = nrb;            /* number rows in b    */
  75.     nfwab = 1;        /* column separation    */
  76.     incrb = ncb;        /* column elem incrmnt    */
  77.     }
  78.     else            /* m1 = ABT or ATBT    */
  79.     {                /* set for b-transpose    */
  80.     mcb   = nrb;        /* n rows -> n cols    */
  81.     mrb   = ncb;        /* n cols -> n rows    */
  82.     nfwab = ncb;        /* column separation    */
  83.     incrb = 1;        /* column elem incrmnt    */
  84.     }
  85.  
  86.     ij = ic = 0;
  87.  
  88.     for (i = 1; i <= mra; ++i)
  89.     {
  90.     ik = 0;
  91.  
  92.     for (j = 1; j <= mcb; ++j)
  93.     {
  94.         ia = ij;
  95.         ib = ik;
  96.  
  97.         /* ----------------------------- */
  98.         /* calculate next vector product */
  99.         /* ----------------------------- */
  100.         r = 0.0;
  101.  
  102.         for (k = 1; k <= mrb; ++k)
  103.         {
  104.         r += a[ia] * b[ib];
  105.  
  106.         ia += incra;
  107.         ib += incrb;
  108.         }
  109.  
  110.         /* --------------------------------- */
  111.         /* Check if negative product desired */
  112.         /* --------------------------------- */
  113.         if (m2 == MP || m2 == CMP)
  114.         r = -r;
  115.  
  116.         /* ---------------------------------- */
  117.         /* Check if replace-add (or subtract) */
  118.         /* ---------------------------------- */
  119.         if ((m2 == P) || (m2 == MP))
  120.         c[ic] = r;    /* Just store product */
  121.         else
  122.         c[ic] += r;    /* replace-add (or subtract) */
  123.  
  124.         ik += nfwab;
  125.         ic++;
  126.     }
  127.  
  128.     ij += nfwaa;
  129.     }
  130. }
  131.  
  132.  
  133.